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Abstract 

A symmetrized Density Matrix Renormalization Group procedure together with the correction vector approach 
is shown to be highly accurate for obtaining dynamic linear and third order polarizabilities of one-dimensional 
Hubbard and U — V models. The U — V model is seen to show characteristically different third harmonic generation 
response in the CDW and SDW phases. This can be rationalized from the excitation spectrum of the systems. 

Keywords: (Semiempirical models and model calculations, Many-body and quasiparticle theories, Excitation spectra 
calculations, Models of non- linear phenomena, Non-linear optical method) 



In the electronic structure theory of interacting quasi- 
one-dimensional systems, the density matrix renormal- 
ization group method (DMRG) has established itself as 
the method of choice. As introduced by White]!], [^], 
this is a static method for obtaining ground and low- 
lying excited states of model systems. We recently ex- 
tended the method to excited states by exploiting sym- 
metries of model systems. In this paper, we further 
extend the method for obtaining accurate dynamical 
linear and nonlinear optic(NLO) coefficients. 

The model exact dynamical NLO properties^, [|, [5| of 
the Pariser-Parr-Pople (PPP) and Hubbard models for 
finite chains were obtained for the first time by Soos and 
Ramasesha by solving for correction vectors The 
two correction vectors (f>[ (u>i) and (jy?) (u>i, u^) encoun- 
tered in the computation of polarizability and third- 
order polarizability are defined by the following linear 
equations. 



can be written as, 
a tJ (cu) =< ^ (1) M|A,|G > + < ^\-uj)\HG >, (3) 

7ii«( w l> W 2,W 3 ) = P < (j)f\(jJa)\Jij\4>^i {-0J1-0J2, -U)i) > 

(4) 

where the operator P generates all the permutations: 
(— u) a ,i), (wi, j), (a>2, k) and (u)3,l) leading to 24 terms 
for ^jijki with u)„ = —lui — L02 — W3. The tumbling aver- 
aged a and 7 can be defined as 

I 3 1 3 



i=l 



{H-E + w 1 +ie)\$\u 1 )=ji i \G> 



(1) 



(H-£ +^2 + *e)|<4 2 Vi,^ 2 ) =Mil4 1} M > ( 2 ) 
where /2,s are the dipole displacement matrices, \G > 
and Eq are the ground state eigenfunction and energy, 
u> is the excitation frequency, e is the damping factor 
and H is the Hamiltonian matrix. The correction vector 
is solved for in the basis of the configuration functions 
(Soos and Ramasesha employed a valence bond basis), 
which is also the basis in which the Hamiltonian ma- 
trix is set-up for obtaining the ground state. In terms 
of these correction vectors, the components of the po- 
larizabilities, «y, and third-order polarizabilities, Jijki, 



to allow comparison of the calculated NLO response 
with experiments on systems containing molecules in 
random orientations Rj. 

The above method is quite general and has been 
employed in the computations of dynamic NLO coef- 
ficients of a wide variety of Hamiltonians B. The in- 
homogeneous linear algebraic equations encountered in 
this method involve large sparse matrices and an it- 
erative small matrix algorithm, which runs parallel to 
the Davidson algorithm for eigenvalue problems, gives 
rapid convergence for the solution of the system of 
equations M. 

The DMRG method [FjJ as implemented, readily gives 
the ground state energy and eigenfunctions and the 
Hamiltonian matrix. The extension of the DMRG 
method to calculate NLO coefficients by the correc- 
tion vector route requires, in addition, the dipole dis- 
placement matrix to be constructed at each iteration. 



Table 1 : Comparison of DMRG and exact a and 7 values at u> = O.lt, e = O.OOli for Hubbard and U — V models 
of 12 sites, a values are in units of 10 -24 esu while 7 values are in units of 10 -36 esu. The all — trans polyene 
geometry is assumed in the calculations and bond length is taken as 1A for the Uniform chain. 



NLO 


U = M, V = 


U = 3t, V = 1.5* 


U = 3t, 


V = 2.0t 


Coefficients 


Exact 


DMRG 


Exact 


DMRG 


Exact 


DMRG 


&XX 


42.61 


42.57 


152.12 


152.05 


997.90 


997.72 


C*xy 


2.95 


2.91 


35.45 


35.38 


383.11 


382.99 


a yy 


1.52 


1.44 


14.20 


14.06 


154.73 


154.41 


a 


14.71 


14.67 


55.44 


55.37 


384.21 


384.04 


'Jxxxx 


26571 


26566 


51380 


51374 


-32929131 


-32929014 


Ixxyy 


113.8 


109.1 


2637.7 


2628.6 


-4915213 


-4915102 


lyyyy 


1.24 


1.22 


189.9 


186.0 


-736663 


-736556 


7 


5367 


5356 


11357 


11345 


-8704799 


-8704671 



The matrices of the dipole operators for the fragments 
are constructed in the DMRG scheme by renormaliz- 
ing the matrix representations of the dipole operator 
corresponding to the left and right parts of the sys- 
tem using the density matrix eigenvector basis in a way 
which is completely analogous to the construction of 
the corresponding Hamiltonian operators for the frag- 
ments. The matrix representation of the dipole opera- 
tors for the full system are obtained as direct products 
of the fragment matrices analogous to the way by which 
the full Hamiltonian matrix is constructed. The dipole 
displacement matrices are obtained by subtracting the 
corresponding components of the ground state dipole 
moment from the diagonal elements of the dipole ma- 
trices. With the computation of the dipole displace- 
ment matrix, the equations for the correction vector 
are completely defined and obtaining the NLO coeffi- 
cients appears to be straightforward. However, unless 
symmetry of the model Hamiltonian is appropriately 
exploited, these equations cannot be solved for frequen- 
cies of interest, in most cases. 

In the case of model Hamiltonians for polymers, the 
systems usually possess total spin symmetry, reflection 
symmetry about the middle of the chain and in some 
cases electron- hole symmetry. These symmetries ensure 
that a given correction vector spans only a symmetrized 
subspace and not the entire Hilbert space of the Hamil- 
tonian of the given system. The correction vector lies 
in the symmetry subspace which is connected to the 
ground state by the electric dipole operator. For u> 
values corresponding to resonance between the ground 
state and eigenstates in that particluar symmetry sub- 
space, Ihs of eqns(l) and (2) become singular for e = 
and present numerical difficulties even when solving 
them for reasonable nonzero e values. However, if we 
do not exploit the symmetries of the Hamiltonian, we 
encounter singularities in eqns(l) and (2) even for those 
to values corresponding to eigenstates of the Hamilto- 



nian found in other symmetry subspaces which are not 
connected to the ground state by the dipole operator. 
Therefore, numerically it would be impossible to obtain 
the correction vectors at these frequencies and thereby 
the associated response of the system. 

For example, in Hubbard chains at intermediate cor- 
relation strengths, a triplet excited state lies below the 
lowest singlet state in the ionic B subspace [0. The 
states in the ionic B space are connected to the ground 
state via one-photon transitions. The resonances in po- 
larizability are thus expected only at frequencies cor- 
responding to the energy levels in the ionic B space, 
relative to the ground state. However, we cannot solve 
for the correction vector using eqn(l) at an excitation 
energy corresponding to the energy of the lowest triplet 
state. Thus, the technique of correction vectors will not 
be able to give the complete dispersion of the polariz- 
abilities upto the first one-photon resonance, unless in- 
terferences due to spurious intruders such as the triplet 
states are eliminated by suitably block-diagonalizing 
the Hamiltonian matrix. This problem of intruders 
becomes more severe with increasing system size due 
to increasing number of states intruding below the fre- 
quency corresponding to the first "true" resonance. 

In the symmetrized DMRG procedure, we have em- 
ployed Parity (invariance of Hamiltonian for inter- 
changing second quantized operators of up and down 
spins for an even electron system in the M s = sec- 
tor), electron- hole symmetry and inversion symmetry of 
the chain. Parity classifies the spin states as odd(o) or 
even(e) corresponding to odd or even integer value of to- 
tal spin. Electron-hole symmetry bifurcates the states 
into ionic{-) and covalent(+) spaces and exists only at 
half-filling. The inversion symmetry of the chains al- 
lows classifing the states as A g or B u . These basic sym- 
metries are sufficient to exclude spurious resonances in 
the solution of eqn(l) and (2) for excitation frequencies 
in the range of interest ( usually u> < E g /2). The im- 



plementation of symmetry within the DMRG procedure 
has been detailed elsewhere [fTl) 

The symmetry adaptation scheme for excited states 
and NLO properties has been implemented both within 
the infinite chain DMRG and finite system DMRG al- 
gorithms. In the latter, we incorporate the spatial sym- 
metry only when the left and right parts of the system 
are identical, i. e. at the end of each finite system itera- 
tion. Without iterating over the density matrices of the 
fragments ( i.e. within the infinite system algorithm), 
we find that the energy difference between a chain of 
length N with N+l and N—l electrons is equal to the 
Hubbard parameter U to an accuracy of w 10~ 3 . Af- 
ter three iterations of the finite system algorithm, the 
accuracy improves to w I CP 5 , for a U value of At. 

To illustrate the power of the symmetrized DMRG 
method, we present results of DMRG calculations for 
ground and excited states of uniform Hubbard chains 
at half-filling, for U/t of 4.0 and 6.0, with chain lengths 
of upto 50 sites. We have obtained the lowest en- 
ergy states in all the eight subspaces of the Schrodinger 
group of the system, retaining 100 to 150 eigenvectors 
of the density matrix at each iteration. The extrapo- 
lated ground state energy per site of the infinite chain 
for both values of U/t agrees with the exact values 
of 0.5737331* and 0.4200125* respectively to 5 decimal 
places. Model exact gaps vanish for infinite chain. 

The energy of the one-photon transition from the 
ground state (lying in the e A + subspace) to the low- 
est energy state in the e B~ subspace defines the opti- 
cal gap of the chain. The extrapolated DMRG optical 
gap compares very well with the optical gap from the 
Bethe ansatz solution for the infinite chain. The ex- 
trapolated infinite chain values from a quadratic fit in 
1/n are 1.278* and 2.895*, obtained with a cut-off of 
m=150 for U/t = 4.0 and m=100 for U/t = 6.0 while 
the Bethe- ansatz values are 1.2867* and 2.8926* respec- 
tively n 

In fig.l, we present the lowest spin excitation gaps 
as a function of inverse chain length for two different 
values of U/t of 4.0 and 6.0 for uniform Hubbard chains 
with upto 50 sites. The lowest triplet state lies in the 
covalent A subspace but with an odd parity. The spin 
gap is defined as the energy gap between the ground 
state singlet and the lowest triplet in the °A + subspace. 
The Bethe- ansatz solution yields a vanishing spin gap 
in the thermodynamic limit for uniform Hubbard chain 
at half filling. Second order polynomial fit in 1/n to 
our DMRG data are consistent with the exact results. 

We compare the polarizabilities and third order po- 
larizabilities from the DMRG calculations with the 
model exact results for the Hubbard and the U — V 
models, in table 1. We find that the DMRG values 
agree very well with the model exact results. 

We have computed the dynamic NLO coefficients of 



theU-V model in the SDW (U > V/2) as well as CDW 
(U < V/2) phases with and without bond alternations. 
In fig. 2 we show the dependence of ln(pf) on the number 
of sites in the chain, N (the latter also on a logarithmic 
scale). For small U/t and S = 0.0, the dependence 
exhibits a power-law behaviour (7 ~ L b , b = 3.62). As 
we increase U /* or 5, the chain length for which the 
power law behaviour is seen is reduced. This shows 
that the 7r-coherence length decreases with increasing 
U /t and S. It is also interesting to note that in the SDW 
limit, the dimerized chain (S = 0.09) exhibits a weaker 
THG response than uniform chain(5 = 0.0). A similar 
behaviour is also seen for linear polarizabilities Jl3| . 

In the CDW phase, the optical gap reduces rapidly 
with chain length. At moderate chain length (N > 10) 
it is observed that the optical gap is small enough that 
the first three photon resonance occurs below fioj = 0.1* 
(Optical gap E g becomes < 0.3*) and this is reflected 
in a change of sign in all components of 7. In the limit 
of very long chains, the optical gap vanishes and we 
should expect 7 to be negative for any finite but small 
excitation frequencies. The behaviour of 7 in the CDW 
phase is qualitatively different from what is observed in 
the SDW phase. In the latter, the optical gap remains 
finite and the 7 is expected to remain positive. 

The static 7 values cannot be computed from correc- 
tion vector method as the equations for (fy^ (u>x , 0J2) be- 
come singular for 0J2 = 0. However, following Soos and 
RamaseshaQ, it is possible to obtain static 7 values 
from the dynamic values at two different frequencies 
by a Taylor expansion. We have computed j X xxx(0) 
using this method from THG coefficients at two dif- 
ferent frequencies of 0.1* and 0.2*, at the crossover re- 
gion (U = 2V) between SDW and CDW phases. In 
fig. 3 we present In ^xxxxiO) vs In L behaviour. The 
Jxxxx(0) value is positive as the thermodynamic limit 
of the CDW phase is not reached even at the maximum 
chain length (20 sites chain) we have studied, jxxxxi®) 
has a power law dependence on chain length L with an 
exponent of 3.94. This is very close to the value seen 
for PPP models in an earlier exact study of finite chains 
in static limit 

To conclude, the symmetrized DMRG method has 
been successfully extended to compute dynamic NLO 
coefficients. This has been illustrated by comparison 
with exact calculation on Hubbard and U — V models. 
The latter has been studied in both the CDW and SDW 
regimes and the qualitatively different behaviour of 7 
in these two regimes is highlighted. 
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Figure Captions 

Fig.l 

Spin Gap(dcfined in the text) as a function of inverse 
chain length 1/N for uniform Hubbard chains with U/t 
4.0 and 6.0. Circles correspond to a DMRG cut-off of 



m = 100 and dots to to = 70. 
Fig.2 

Log-log plot of 7 with site number. The model param- 
eters are (i)U = 3t,V — i(squares), (ii)U = it, V = 
t,S = 0.09(trianglcs), (iii)U = At,V = t(diamonds) and 
(iv)t/ = At, V = t, S = 0.09(circles). 
Fig.3 

Log-Log plot of 'fxxxxi®) with chain length for the pa- 
rameter U = 3t, V = 1M, S = 0.09. 



